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The hyperfine structure (hfs) and the g factor of a bound electron are caused by external magnetic 
fields. For the hfs, the magnetic field is due to the nuclear spin. A uniform-in-space and constant-in- 
time magnetic field is used to probe the bound-electron g factor. The self-energy corrections to these 
effects are more difficult to evaluate than those to the Lamb shift. Here, we describe a numerical 
approach for both effects in the notoriously problematic regime of hydrogen-like bound systems 
with low nuclear charge numbers. The calculation is nonperturbative in the binding Coulomb field. 
Accurate numerical values for the remainder functions are provided for 2P states and for nS states 
with n = 1, 2, 3. 
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I. INTRODUCTION 

The interaction of a bound electron and an atomic nu- 
cleus is characterized by the parameter Za, where Z is 
the nuclear charge number and a is the fine-structure 
constant. This universal "coupling parameter" sets the 
scale for calculations of the radiative corrections to var- 
ious bound-state effects including the hyperfine struc- 
ture (hfs) and the bound-electron g factor. Traditionally, 
theoretical investigations of radiative corrections in light 
systems relied upon an expansion in powers of Za and 
\n(Za). However, today it is desirable to advance theory 
beyond the predictive limits given by the highest avail- 
able terms in the Za-expansion. This can be done by 
carrying out calculations with using nonperturbative (in 
Za) propagators. Such calculations demand rather so- 
phisticated numerical techniques, which were developed 
relatively recently. Indeed, all-order calculations of the 
self-energy (SE) correction in the presence of a magnetic 
field started in the 1990s [E 0, i, 11, H, 0] , extending dur- 
ing past years to a wide range of reference states and 
nuclear charge numbers @, i, [|, El E2, El El El • 

Numerical calculations of the SE corrections are par- 
ticularly difficult for low values of Z. This is mainly 
for two reasons. First, the goal of the calculations is 
the contribution beyond the known Za-expansion terms. 
For the hfs, the higher-order effects are suppressed with 
respect to the leading correction by a factor of (Za) 3 . 
For the g factor, they enter only at order of (Za) 5 and 
thus become very small numerically in the low-Z region. 
Second, in actual calculations there are additional can- 
cellations arising at intermediate stages of the numerical 
procedure. These cancellations become more severe for 
smaller values of Z and lead to further losses of accuracy. 

In this article, we treat the two most important exam- 



ple cases of the bound-electon SE corrections in external 
magnetic fields: the SE correction to the hfs and the 
SE correction to the bound-electron g factor. We evalu- 
ate both of these corrections for the ground and excited 
states of hydrogen and of light hydrogen- like ions. The 
first attempt at an all-order evaluation of the SE correc- 
tion to the hfs of hydrogen was made in Ref. [J] . Because 
of insufficient numerical accuracy, the goal was reached in 
an indirect way: the known terms of the Za expansion 
were subtracted from numerically determined all-order 
results for Z > 5, and the higher-order remainder was 
extrapolated down toward the desired value Z = 1. The 
accuracy of the numerical evaluation of the SE correction 
to the hfs was improved by several orders of magnitude 
during the past years @, E2|- However, the precision ob- 
tained was still insufficient for a direct determination of 
the higher-order SE remainder at Z = 1, and an extrap- 
olation procedure had to be employed again. 

The studies I El reported results for the higher-order 
contribution for the normalized difference of the 15 and 
25 hfs intervals in 3 He + and demonstrated a 2a deviation 
of the theoretical prediction from the experimental result 
[TEl EB[- The accuracy of the extrapolation procedure of 
RefsTyl, E3|, however, has recently become a subject of 
some concern. In particular, an opinion was expressed in 
Ref. [TtJ that the uncertainty of the extrapolation pro- 
cedure should have been estimated as four times larger 
than given in Refs. d El; which would have brought 
theory and experiment into agreement. 

In our recent investigation [lc| , we performed the first 
direct, high-precision theoretical determination of the 
higher-order remainder of the SE correction to the hfs 
of 15 and 25 states of hydrogen and light hydrogen-like 
ions. Good agreement was observed with the previous ex- 
trapolated values [1 E3 > but the accuracy was increased 
by several orders of magnitude. In the present paper, we 
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report the details of this calculation and extend it to the 
higher excited states (35, 2Pi/ 2 , and 2P 3 / 2 ). 

The SE correction to the bound-electron g factor is 
of particular importance because it is used in the de- 
termination of the electron mass value from the experi- 
mental results for the g factor of light hydrogen-like ions 
[l9| . Already at the present level of experimental accu- 
racy, calculations of the bound-electron g factor should 
be performed to all orders in Za. A number of all-order 
evaluations of the SE correction to the a fa ctor have been 
accomplished during last years 0, H, l7Lll(A [ITj | , which re- 
sulted in an improvement of the precision of the electron 
mass value. However, in order to match the 10 -12 level 
of accuracy anticipated in future experiments on the he- 
lium ion [2fJ, the precision of numerical calculations of 
the SE correction should be enhanced by several orders 
of magnitude. 

First results of our evaluation of the SE correction 
to the bound-electron g factor for the 15* state of light 
hydrogen-like ions were reported in Ref. [Hj]. In the 
present investigation we extend our calculation to the 
higher excited states (25, 35, 2P 1 / 2 , and 2P 3 / 2 ) and to a 
wider region of the nuclear charge number Z. Relativis- 
tic units (h = c = m = 1) and Heaviside charge units 
(a = e 2 /47r, e < 0) are used throughout the paper. 

Our investigations are organized as follows. In Sec. [ill 
we discuss general formulas pertaining to the formulation 
of the effect within the formalism of quantum electrody- 
namics. We continue with a detailed description of the 
numerical approach in Sec. IIIIl Numerical results are 
presented in Sec. IIVI We conclude with a summary in 
Sec. El 
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FIG. 1: Feynman diagrams representing the SE correction 
in the presence of an external perturbing field. The double 
line indicates the bound electron propagator, which is non- 
perturbative in the coupling parameter Za and entails an 
arbitrary number of Coulomb interactions with the atomic 
nucleus. The wavy line that ends with a cross denotes the 
interaction with the perturbing potential SV. The latter is 
given by the magnetic field of the nucleus in the case of the 
hfs and by a constant external magnetic field in the case of 
the bound-electron g factor. 



II. GENERAL FORMULAS 

The SE correction in the presence of a binding 
Coulomb field and an additional perturbing potential 
SV is graphically represented by the Feynman diagrams 
shown in Fig. [TJ The general expression for them can be 



conveniently split into three parts [2l| . 

A£ S E = A^ir + A^rod + A£ v . 



(1) 



which are referred to as the irreducible, the reducible, and 
the vertex contribution, respectively. 

The vertex contribution is induced by the diagram in 
Fig.QTb). It can be expressed as 



A-Ever — 7T~ 
Z7T 



E 



du 



(7111(51/1712) {an 2 \I(w)\nia) 



u (l-iO)][e„-w 



,(1-20)] 



(2) 



Here, / is the operator of the electron-electron interaction 

J(w) = e 2 a^a„ (uj) , (3) 

where D^ v is the photon propagator and a M = (1, a) are 
the Dirac matrices. The sums over n\ and 712 involve 
both the positive-energy discrete and continuous spectra 
and the negative-energy continuous spectrum. 

The irreducible contribution is induced by a part of 
the diagrams in Fig. Ufa) and (c) that can be expressed 
in terms of the first-order perturbation of the reference- 
state wave function by SV, 



\6a)= £ 

n 



i)(n\6V\c 



(4) 



The expression for the irreducible contribution is 

AE ir = (H 7 °E( £a )| a ) + (a|7°£(£a)|<5a) , (5) 

where X = S — 5m, 8m is the one-loop mass counterterm, 
and S is the one-loop SE operator, 

/>oo 

T,(e,xi,x 2 ) — 2iaj° / dtoa^ 



xG(e - u, xi,x 2 ) a v D^(u, x 12 ) .(6) 

In the above, G denotes the Dirac Coulomb Green func- 
tion G(e) = [e — H(l — z0)] _1 , H is the Dirac Coulomb 
Hamiltonian, and X\ 2 = x± — x 2 . 

The reducible contribution is induced by a part of dia- 
grams in Fig.QJa) and (c) that can be expressed in terms 
of the first-order perturbation of the reference-state en- 
ergy. It reads 



AE rcd = 5e a (o|t° £s(e) 
ae 



(7) 



where Se a — (a\SV\a). 

Up to now we did not specify the particular form of the 
perturbing potential SV, assuming only its locality. In 
the following, we consider two particular choices of SV, 
both representing interactions with the magnetic field: 
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the Ms interaction and the interaction with the magnetic 
field (the Zeeman effect). In the case of the hfs interac- 
tion, the perturbing potential has the form of the Fermi- 
Breit interaction 



Vfb(t-) 



\e\ a ■ [fj, x r] 

47r r 3 



(8) 



(where fj, denotes the nuclear magnetic moment) and the 
reference-state wave function |a) is the wave function of 
the coupled system (electron-l- nucleus) , 



In the following, all explicit formulas for individual 
contributions will be presented for the case of the hfs. 
When working in the coordinate representation, the cor- 
responding formulas for the g factor can be obtained by 
an obvious substitution. In momentum space, the for- 
mulas for the hfs and for the g factor are different. Our 
present approach to the evaluation of the SE correction 
to the g factor closely follows the one of Ref. [ll| and is 
therefore not described separately. 



\FMpIj) 



E 

Mim a 



C 



FM f 
IMn 



\JMj)\j a m a ), (9) 



where \IMi) denotes the nuclear wave function, \ j a m a } is 
the electron wave function, F is the total momentum of 
the atom, and Mp is its projection. The nuclear variables 
can be separated out by using the standard technique of 
the Racah algebra. It can be demonstrated [l2| that the 
general formulas {J)-® yield contributions to the hfs 
if one employs an electronic perturbing potential of the 
form 



SVhUr) 



E F [r x a] z 



Cms 



(10) 



and takes the reference-state wave function to be the 

electronic wave function with the momentum projection 
l 

2 ■ 



/?? 



(11) 



In the above, Ep denotes the nonrelativistic limit of the 
expectation value of the Fermi-Breit operator on the ref- 
erence state and the prefactor Cms is given by 



C his = m 2 {Zaf 



sign(K J 



(12) 



where n a is the Dirac quantum number of the reference 
state and n a is its principal quantum number. 

In the case of the Zeeman splitting, the perturbing 
potential is 



V Zee (r) 



-e a 



Mr) 



(13) 



where A(r) = i [B X r] is the vector potential. In prac- 
tical calculations, corrections to the Zeeman splitting are 
convenienly expressed in terms of corrections to the g fac- 
tor. It can be easily shown (see Ref. [Tl| for details) that 
the general formulas CO)-® yield contributions to the 
electronic g factor if one employs a perturbing potential 
of the form 



SV g (r) = 2m [r x a]. 



(14) 



and the reference-state wave function with the momen- 
tum projection m a — 1/2. The ^-factor perturbing po- 
tential (TTJ| differs from the hfs potential (|10p only by the 
power of r and the prefactor. 



III. DETAILED ANALYSIS 
A. Orientation 

The general formulas presented in the previous section 
for individual contributions are both ultraviolet (UV) 
and infrared (IR) divergent. In order to obtain ex- 
pressions suitable for numerical evaluation, a careful re- 
arrangement of contributions is needed, together with a 
covariant regularization of divergences. The calculation 
of the irreducible contribution ([5]) can be reduced to an 
evaluation of a non-diagonal matrix element of the first- 
order SE operator ([B]). Its rcnormalization is well known 
and does not need to be discussed here. The numerical 
evaluation of the irreducible contribution was performed 
by a generalization of the approach of Refs. [U, [H[ , with 
the use of a closed-form analytic representation of the 
perturbed wave function \Sa) obtained in Ref. [24[ (see 
also Ref. [U). 

The evaluation of the reducible and the vertex con- 
tribution is carried out after splitting them into several 
parts, 



AE Icd = AE 



(a) 
red 



AE, 



(0) 
red 



AE, 



(i+) 

red 



AE VCI = AE&} + AM°) 



i — ^-^vcr 



(15) 

Ai4? r +) , (16) 



where the upper index (a) labels the contributions in- 
duced by the reference-state part of the electron propa- 
gators and the other indices specify the total number of 
interactions with the binding field in the electron propa- 
gators [the index labels the terms generated by >i 
such interactions!. 



B. Reference-state contribution 



The reference-state contributions AE^\ and AEitl are 
separately IR divergent. The divergences disappear when 
the contributions are regularized in the same way and 
evaluated together. Let us now demonstrate the cancel- 
lation of the IR divergences and obtain the finite resid- 
ual. The part of the vertex and reducible contributions 
induced by the intermediate states degenerate in energy 
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with the reference state is 



1 



(w - iO) 2 



[ ^ (o'|*V|a")(aa // |J(w)|o / o) 

WW 

^(a|<$y|a)(aa'|J(w)|a'a) 



(17) 



where a is the "true" reference state and a' and a" la- 
bel the intermediate states that are degenerate with the 
reference state in energy and have momentum projec- 
tions fi a i and n a ii , respectively. (The intermediate states 
degenerate with the reference state in energy but of op- 
posite parity do not induce any IR divergences because 
of the orthogonality of the wave functions. However, in 
practical calculations we find it convenient to treat all 
the degenerate states on the same footing.) For simplic- 
ity, we now consider the photon propagator in the Feyn- 
man gauge. Then, the operator of the electron-electron 
interaction I takes the form 



where 

D(u,xi 2 ) 



= a a^a^ D(lj, xu) , 

exp(ifc • X12) 



-Am 



dk 



(2tt) 3 w 2 



k 2 



(18) 
(19) 



with p being the photon mass, which regularizes the 
IR divergences. It can be seen that all divergences in 
Eq. (fTT|) originate from an integral of the form 



(w - zO) 2 



(20) 



We now substute Eq. (fl9|) into the above expression, 
twice perform an integration by parts, evaluate the to 
integral by Cauchy's theorem, and obtain 



J = 

7T 



dk _ 



COS kxi2 



(21) 



Adding and subtracting cos k in the numerator of the in- 
tegrand, we separate the above expression into two parts, 
the first of which is convergent when p, — > while the 
other (divergent part) does not depend on x 12 - Setting 
p = in the convergent part and evaluating the integral, 
we obtain 



J = — In £12 

7T 7T 



dk- 



cos k 



Vk 2 



(22) 



The divergent part of J does not depend on the radial 
variables and, being substituted into Eq. ([17]), leads to a 
vanishing contribution. We thus obtain 



AE^ =-[Yl {a'\8V\a"){aa"\a^ \nx l2 \a'a) 



^^{a\8V\a){aa'\a^a^ lna:i2|a'a) 



(23) 



We note that in the case when the perturbing potential 
6V is spherically symmetric, the reference-state contri- 



bution AE^ vanishes as (a'\SV\a") 



J P a "Pa 



(a\SV\a) 



In our case, however, SV represents an interaction with 
the magnetic field, so that AE^ induces a finite contri- 
bution. 

In our practical calculations, the reference-state contri- 
bution was separated from the vertex and reducible parts 
by introducing point-by-point subtractions from the elec- 
tron propagators in the integrands and was calculated 
separately according to Eq. ((23]) . 



C. Zero-potential parts 

The zero-potential parts AE 1 ^ and AEicl are sep- 
arately UV divergent. They are covariantly regular- 
ized by working in an extended number of dimensions 
[D = 4 — 2e) and calculated in momentum space. The 
elimination of UV divergences in the sum of the reducible 
and the vertex contributions is well documented in the 
literature (see, e.g., Ref. so here we operate with 

the renormalized SE and vertex operators, assuming that 
all UV divergences are already cancelled out. 

The zero-potential contribution to the reducible part 
is simple. It is given by 



AE, 



(0) 
reel 



(a\6V\a) 



MP), (24) 



where tp = tp^j is the Dirac adjoint. The derivative of 
the renormalized free SE operator T,^ can be expressed 
as a linear combination of 3 matrix structures, j> = 7^f>^, 
7 , and the unity matrix /, 



dp 



a 

47T 



«i (p) + 7o a 2 (p) + Ia 3 (p) 



(25) 



where p = (m 2 — p 2 ) /m 2 = (m 2 — e 2 + p 2 ) /m 2 and dj (p) 
are scalar functions, whose explicit expression is given 
by Eqs. (53)-(55) of Ref. [2q |. Integrating over angular 
variables, we immediately have 



AE, 



2 = ww> (-£) 



Pr d Pr 

(2tt) 3 



X {ai(p)[£a{g 2 a + fa) + Zprgafa] 

+ a2(p)(g 2 a + f a ) + a 3 (p)(g 2 a - fl)) , (26) 

where p r = \p\ and g a = g a {p r ) and f a = f a (p r ) are the 
upper and the lower components of the reference-state 
wave function in the momentum space. 

The zero-potential vertex part of the SE hfs correction 
is induced by the hfs potential SVm s inserted in the free 
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SE loop. The Ms potential ([10]) in the momentum space 
takes the form 



c-rr f \ E F [q x a] z 

SVhSsiq) = -ft- (-4tt«) ■ 



Cms " q z 
The zero-potential vertex part is then given by 

A £(°) = El (-47TZ) / Eel. Ell 

vcr C Ms [ 4m) J (2tt)3 ( 27 r)3 



(27) 



[q x ^(^1,^2) 



■MP2), (28) 



where q = pi ~p 2 , Pi and j>2 are 4- vectors with the fixed 
time component pi — (s a ,Pi), P2 — (e a ,P2), and Tji is 
the renormalized one-loop vertex operator. For evaluat- 
ing the integrals over the angular variables in Eq. (|2"5)) . 
it is convenient to employ the following representation of 
the vertex operator sandwiched between the Dirac wave 
functions 



i> a {Pi)r R (pi,p2)ipb(P2) = j- 

X [^lX«„A.a(P 1 )°'X-«.M«(P2) 
+ ^2xl KaAl(l (Pl)o-X« aMa (P2) 

+ [KaPl +^4P2)X« /*a(Pl)X«a|*a(P2) 

+ (^ 5 pi +^ 6 p 2 )X-„ a p a (pl)X- KaAta (P2) 



(29) 



where the scalar functions 7?.i = lZi(pi r ,P2r, q_r) are given 
by Eqs. (A7)— (A12) of Ref. [H, and Pl = 
Pir = \Pi\i and g r = \q\. The dependence of the inte- 
grand of Eq. (j2"5|) on the angular variables can now be 
parameterized in terms of the basic angular integrals K{ 
introduced and evaluated in Appendix The result is 



ae\?1 = 



n « s I dpi r dp2r / dg r 

C h fs 48tt 5 J J\ Plr - P2r \ Qr 



x |[-pi r Xi(/c a ) + p 2r K[(K a )]K 1 

+ [-p lr Ki(-K a ) +P2rK' 1 (-K a )}K 2 
~ PlrP2rK 2 (K a ) + T^i) 

lrP2rK 2 (-K a ) (TZ 5 +7?. 6 )| • 



Plr 



(30) 



The above equation was used for the numerical evalua- 
tion. It contains four integrations (the fourth one, over 
the Feynman parameter, is implicit in the definition of 
the functions Ri). All the integration were performed us- 
ing Gauss-Legendre quadratures, after appropriate sub- 
stitutions in the integration variables. We note that the 
integration variables pi r , p 2r , and q r resemble the well- 
known perimetric coordinates (27| , in the sence that they 
weaken the (integrable) Coulomb singularity of the inte- 
grand at q r — 0. 



D. One-potential vertex part 



The one-potential hfs vertex part AEl]) T is given by 



{p" x A(p,pV% 



(q — p") 2 p 



II \2 „/' 2 



■Mp') 



(31) 



where q = p — p' is the total momentum transfer (final 
minus initial) for the electron vertex function A, and the 
time component of the 4-vectors is fixed by Po — p'q — e a 
and p'q — 0. The 4-point vertex function A is given by 

, , „, 16tt 2 f d 4 k 1 
A j (p,p>, p ") = — I 7 —^ [ - (32) 



i J (2tt) 4 k 2 
+ w)7o(^ + m)7 3 -(^ / - jH- m) 7 ' T 

[(p — k) 2 — m 2 ] [{p — k — p") 2 — m 2 ] \{p' — k) 2 — m 2 ] 

The evaluation of AEill is performed by using the stan- 
dard technique for the evaluation of Feynman diagrams 
(for a short summary of the relevant formulas, see Ap- 
pendix D of Ref. [28|). First, we use three Feynman pa- 
rameters in order to join the 4 factors in the denominator 
of the integrand in Eq. (|3"2")) . Denoting the numerator as 
Nj(k), we obtain 



Aj(p,p\p") 



dx dy dz 6x 2 y 
16tt 2 f d*k 



N 3 (k) 



(2tt) 4 [{k - xb) 2 - xA] 4 



(33) 



where x, y, and z are the Feynman parameters (here and 
below it is assumed that all integrals over the Feynman 
parameters extend from to 1). We denote b = (1 — 
y)p + yp' + yzp", and A = xb 2 + m 2 — (1 — y)p 2 — yp' 2 — 
yz(p" 2 + 2p'-p"). 

Next, we shift the integration variable k —>■ k + xb and 
perform the integration over k. The result is 



kj(p,p',p") = dxdydzy 



Njjxb) xN^g^ 



A 2 



2A 



(34) 



where N 2 j is defined so that N^jk^k^ is the quadratic 
in k part of Nj(k). We note that, after shifting the inte- 
gration variable, only even powers of k yield a non-zero 
contribution to the integral (i.e., the terms proportional 
fc M and kpkykp vanish). 

Next, the integration over p" is carried out. We intro- 
duce the function by 

p'l K^p' lP ") 



dx dy dz y 



3 p" 2 (q-p") 2 
dp" 



1 



A 2 



- 2x 



(2tt) 3 p" 2 (q-p") 2 
P'lN^ip") 



A 



(35) 
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where N od (p") = Nj{xb) and N 2 ^(p") = N^g^/4. The 
integral over p" in Eq. ([3!>[) can be expressed in terms of 
the Lewis integral [29j |. However, we prefer to perform 
this integration straightforwardly by merging denomina- 
tors using Feynman parametrization. In this way, we end 
up with an additional integration to be performed nu- 
merically, but the structure of the expressions involved 
becomes somewhat simpler. 

Let us illustrate the further evaluation by considering 
the contribution induced by the first term in the square 
brackets in Eq. ([3"5)) . which will be denoted by So,y . We 
merge the denominators by introducing two more Feyn- 
man parameters, 



where w = 1 — xyz, 



du dt 



6u 2 t 



(wyz) 2 \{p" - uc) 2 + ufl} 4 ' 
(36) 



c = -[x(l - y)p - (1 - xy)p'} + (1 - t)q , (37) 
w 



and 



Q = -uc 2 + (1 - t)q 2 + 



wyz 



x - y)p + yp'f + m 2 - (1 - y)p 2 - yp' 2 } . 

(38) 

Substituting Eq. ([3u]) into Eq. ([33)) and shifting the inte- 
gration variable, we get 



So,iiO,p') = J di 



6u 2 t 



dp" 



Mi, 



, (39) 



yw 2 z 2 J (2tt) 3 (p» 2 + U Q) 4 

where dp = dx dy dz du dt and 

M i:j = [p'l + uci) N 0J (p" + uc) 
ee M 0)ij + M^p'l + Mf^p'ip'l 

+ M^plp'lpl % + Mf^plp'lp'^pl . (40) 

The above equation defines the M functions as the 
coefficients from the expanded form of the expression 
{p'{ + uci ) Nqj (p" + uc) . Performing the integration over 
p" in Eq. ([311), we obtain 



3o,«(p,j/) =g^/fc 



u 2 t 



I. 



3M 



0,ij 



yw 2 z 2 \ (ufl) 5 / 2 



(ur>) 3 / 2 



{ u ny/ 2 



(41) 



Because VI is linear in u, the integral over u is elementary 
and can be expressed in terms of logarithms. The four 
remaining integrations over the Feynman parameters re- 
main to be evaluated numerically. To complete the eval- 
uation of Sq^j, one needs to obtain explicit expressions 



for the numerators M/^- and to bring them to the stan- 
dard form. Under "the standard form" we understand a 
linear combination of independent matrix structures, see 
below. This is the most tedious part of the calculation 
since the expressions involved are very lenghty. Usage of 
symbolic computation packages is indispensable in this 
case. 

Having obtained an expression for Sy, we write the 
correction to the hfs as 

x S(p r ,j£,? r ;* 1 ,... J Xaa)Va(p / ), ( 42 ) 

where we used the notation 5 = eoij Hy with ey/- de- 
noting the Levi-Civita symbol. In Eq. (|42[) . we indicate 
explicitly the dependence of S on 32 basic matrix struc- 
tures Xi. The main four of these are: [p x -f] z , [p' x j] z , 
[pxp'] z , and [7 x 7] z . The rest is obtained by multiplying 
each of them by ^, p 7 , j>f , 7 , ^7°, and ^7°^'. 

In order to perform the integration over all angular 
variables in Eq. (|42[1 except for £ = p ■ p' , we define the 
angular integrals Yi that correspond to the basic matrices 
Xi by 



dp dp' ip a (p) XiF(p r ,p' r ,q r )ip a (p') 



(43) 



J d£YiF(p r ,p' r ,q r ) , 



?: = !,..., 32. 



where F is an arbitrary function. All Yi may be expressed 
in terms of the elementary angular integrals listed in Ap- 
pendix rAj 

Using the angular integrals Yi, we can write the final 
expression for the one-potential vertex term suitable for 
a numerical evaluation, 



x / dq r prp' r qrZ(p r ,p' r ,q r ;Y 1 ,...,Y 32 ). 

J |Pr— Pr\ 

(44) 

Altogether, Eq. (j44|) contains 7 integrations to be per- 
formed numerically, 3 of them being written explicitly 
and 4 Feynman-parameter integrations contained in the 
definition of the function S. The numerical evaluation 
was performed using Gauss-Legendre quadratures for all 
integrations. In order to prevent losses of accuracy due 
to numerical cancellations, we used quadruple-precision 
arithmetic (accurate to roughly 32 decimals) in a small 
part of the code, which was identified to be numerically 
unstable. The evaluation was rather time-consuming 
(about a month of processor time for each value of Z 
and each state) and was performed with the help of the 
parallel computational environment at MPI Heidelberg. 

The one-potential vertex part has been crucial to our 
calculation, and so it may be appropriate to summarize 
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once more the basic steps in its evaluation: First of all, let 
us recall that our "one-potential vertex part" actually in- 
volves two vertices inside the loop, one being a Coulomb 
vertex and the other being a magnetic vertex (coupling 
to the external field). Therefore, there are three fermion 
propagators inside the loop and one photon propagator, 
necessitating the introduction of three Feynman param- 
eters to join denominators. The incoming Coulomb mo- 
mentum and the exchanged momentum with the exter- 
nal field entail two further Feynman parameters, one of 
which is integrated out analytically. In addition to the 
four remaining Feynman parameters, we have two radial 
integrations over the absolute values of the initial (p 1 ) 
and final (p) electron momenta, and an integration over 
the direction cosine £ (transformed by a change of vari- 
able to an integration over q r = \p — p'\). The three 
additional integrations account for the resulting seven- 
dimensional integral. In the corresponding calculation in 
free QED, one could hope to carry out the radial inte- 
grations analytically, because the incoming and outgoing 
fermions are on the mass shell and described by plane 
waves. Here, however, the bound states are being off 
the mass shell and have a much more complicated struc- 
ture, so that the radial integrations have to be evaluated 
numerically. The separate calculation of the full one- 
potential vertex part as described in the current section 
leads to a numerically favourable scheme, because this 
part can be then subtracted from the integrand of the 
remaining nonperturbative vertex contribution, thereby 
leading to a drastic improvement in the convergence of 
the resulting partial- wave expansion (see Table U below) . 



E. Many-potential vertex part 

The general expression for the many-potential vertex 
part AEitt is obtained from Eq. by applying the ap- 
propriate set of subtractions in the electron propagators. 
The required subtractions are given by 



GSVG -> G 5V G - G {a) SV G (q) 

_ G (o) 6VG (o)_ G (o) SVG (i) 



(45) 



G (D SVG W 



where G denotes the bound-electron propagator, G^ is 
the free-electron proparator, is the electron prop- 
agator with one interaction with the binding Coulomb 
field, and G^ is the reference-state part of the bound- 
electron propagator. This subtraction takes into account 
all terms which have been calculated separately using 
different approaches, as described above. 

(2+) 

In order to perform a numerical evaluation of A_Ever , 
it is convenient to rotate the integration contour of the 
photon energy u) from (— oo,oo) to be parallel to the 
imaginary axis of the u) complex plane. In this work, 
we define a deformed to integration countour Clh con- 
sisting of two parts, a low-energy part Cl and a high- 
energy part Cjj. The low-energy part contains the in- 
terval we(A- iO, —iO) on the lower bank of the cut of 



the photon propagator and the interval (iO, A + iO) on 
the upper bank of the cut, with A = Zae a . The high- 
energy part consists of two intervals, (A + iO, A + ioo) 
and (A — iO, A — ioo). The contour Clh defined in this 
way differs from the one used by P. J. Mohr [3Cj only by 
the choice of the separation point A (the value A = e a 
instead of A = Zae a was employed in Ref. [30|). 

(24-) 

The high-energy part of A£ver is given by 



AE, 



(2+) 
■ H 



1 

-Re 

7T 



n\n 2 _ 



(ni\SV\n2) (aii2\I(A + iu)\nia) 



{e a - A-ioj- e ni )(e a - A 



iuj — e, 



subtractions 



(46) 



where the subtractions are given by Eq. (|45p . The low- 
energy part of AEva needs a careful treatment because 
of single and double poles situated near the contour Cl, 
which are due to virtual bound states of lower energy 
than the reference state. The single poles can be in- 
tegrated via a Cauchy principal value prescription, and 
the double poles can be converted to single poles via an 
integration by parts. We thus write the low-energy part 

of AMer +) as 



k Jo 



E 



0<£„<e a 



E 

•1 = ' 

F' M 

£ a - tO - £r. 



F ni n 2 (u>) 



111 ™2 

not < s ni — e n2 < e a 



(e a - uj - e ni )(e a — a — e„ 2 ) 



subtractions 



F nn {A) 



K n^~S £a - A - £„ ' 

0<E T ,<£ a 



where 



F niri2 (uj) = (ni\SV\n 2 ) (<m 2 |Im [J(w)] \nia) 



(47) 



(48) 



the prime denotes the derivative over u), and P denotes 
the Cauchy principal value of the integral. In Eq. (|47[) . 
all terms that induce double poles on the interval lu G 
(0, A) (i.e., intermediate states with < e ni = e n . 2 < e a ) 
have been integrated by parts. We recall that the term 
with e ni — £n 2 — £ a is removed by the G^ part of the 
subtraction (|45]) . 

The need to evaluate the principal value of the integral 
over lu complicates the numerical calculation of the low- 
energy part. In the case when there is a single pole only 
(which takes place for the 2s and 2pi/ 2 reference state), 
the problem is most easily solved by employing a numeri- 
cal quadrature symmeric around the position of the pole. 
In the general case with more than one singularity to be 
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treated, this approach is not effective. A better way is to 
introduce subtractions in the integrand that remove the 
singularities at the poles and to evaluate the principal 
value of the integral of the subtracted terms analytically. 
We introduce the subtractions by observing that the fol- 
lowing difference does not have any singularities on the 
interval u> € (0, A), 



E 



not < £ ni — Snr, < Sa 



(m |(5l / |n 2 ) (an 2 |Im [I(to)] |ma) 
(e a - to - s ni )(e a - u — £„ 2 ) 



E 

0<£„ 1 <£ c 



(a(5ni|Im[/(e a - e ni )] \n\a) 



E„ — tO — £r. 



(an 2 |Im[7(e a - e» 2 )] \Sn 2 a) , 



0<£„,<e a 



GO - E T , 



where 



\8m) = J2 



\n 2 )(n 2 \SV\ni 



(50) 



We note that the terms with (e ni = e a ,e n2 ^ e a ) and 
(e ni ^ e a , £ n2 — £ a) present in Eq. (|49|) do not induce any 
singularitites because Im[/(0)] = 0. The perturbed wave 
function \8rii) is known analytically (for the hfs perturb- 
ing potential, both the diagonal and the non-diagonal in 
k parts; for the ^-factor perturbing potential, only the 
diagonal in k part) from the generalized virial relations 
for the Dirac equation [24|, [H[ . 

In order to complete our discussion of the evaluation of 
the many-potential vertex part, we present the explicit 
expression for it after the integration over the angular 
variables. This expression reads 



AE (2+) = E F "« 



Chfs 2-7T 



? / d " T 

L H 





J2 




{ Ja 


ja 


1} 



P(ni, n 2 ) Rl (oo, an,2nia) 



{Ea 



£«i)(£a — tJ — £n 2 ) 



subtractions 



(51) 



where Rl is a relativistic generalization of the Slater 
radial integral, whose explicit expression is given in 
Ref. [321 . P(ni,n 2 ) is given by 

p{m,n 2 ) =(-iy« +1/2 c}°_ 1/2jal/2 

x K1 + K2 (_ K2 || C7 (i)|| Kl ) fl _ 2 („ 1)n2 ) ) (52) 

where C^™ ■ m is the Clebsch-Gordan coefficient, 

Cm = \j 47r/ (21 + 1) Yim is a normalized spherical har- 
monic, and 



R- 2 (rii,n 2 ) = / dr[g ni (r)f ri2 (r) + f ni (r)g n2 (r)] 



(53) 



The numerical evaluation of the many-potential vertex 
contribution is the most difficult part of the calculation. 
The key feature that limits the accuracy achievable in a 
numerical calculation is the convergence of the partial- 
wave expansion. We recall that the many-potential ver- 
tex contribution AE^cr contains two and more Coulomb 
interactions (and a magnetic interaction) inside the self- 
energy loop. The convergence of its partial-wave expan- 
sion is much better than that for the vertex contribution 
with just one Coulomb interaction. In order to illustrate 
this point, Table [I] presents a comparison of the partial- 
wave expansion of the vertex contribution with two and 
more Coulomb interactions, Ai?vcr , and of that with one 
and more Coulomb interactions, AEic^ . It can be seen 
that the subtraction of the one-potential vertex contribu- 
tion improves the numerical accuracy by about 5 orders 
of magnitude. 

The partial-wave expansion was cut off at the maxi- 
mum value of |K m ax| = 120. Quadruple-precision arith- 
metics was required for the Dirac Green function in or- 
der to control the numerical accuracy at the required 
level. This computation was performed with a quadruple- 
precision generalization of the code for the Dirac Green 
function developed in Refs. [1, [Hj]. It should be men- 
tioned that the evaluation of the high-energy part of 

(2+) 

A£vcr with the integration contour Clh requires the 
Dirac Green function with general complex values of the 
energy argument. (This is in contrast to the approach 
used in Refs. [H, [H, where the integration contour 
is chosen in such a way that only the real and purely 
imaginary values of the energy argument are required.) 
The computation of the Dirac Green function for general 
complex energies to becomes numerically unstable when 
k is large and arg(w) is close to 7r/4. Because of this, 
we were not able to extend the partial-wave summation 
beyond |K max | = 120. 

The general scheme of our evaluation is as follows. We 
perform the summation over k directly in the integrand, 
before any integrations. The summation is terminated 
when a suitable convergence criterion is fulfilled or when 
the cutoff value |ft ma x| is reached. In order to estimate 
the dependence of the final result on the cutoff param- 
eter, results for several intermediate cutoffs are stored, 
each consequent one being twice larger than the previ- 
ous (see Table [I] for an illustration). The omitted tail 
of the expansion was estimated by using the e algorithm 
for Pade approximation, and the uncertainty of the ex- 
trapolation was taken about 50%-200% of the estimated 
tail. 



F. Many-potential reducible part 

According to Eq. ([7]), the reducible part of the SE cor- 
rection involves the derivative of the SE operator, "sand- 
wiched" in the reference state. The zero-potential part 
of the reducible contribution has already been discussed 
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TABLE I: Comparison of the convergence of the partial-wave expansion for the corrections AEilt and AE^t for the hfs 
of the 15* state of atomic hydrogen (Z = 1), in units AE/[a/ir Ep]. S(n max ) is the sum of all partial contributions with 
|k| < ftmax, and the convergence is measured as K max is increased. 8S is the increment. For AE^\ at the same value of K max , 
the apparent convergence gives us roughly five more decimals as compared to AE^lf - 1 . 

Af*+) Ap( 2 +) 

L -*- LJ vcr *— i - L/ ver 



^max 


5S 


S ( K max ) 


5S 


S'(^max) 


3 


1.580875 


1.580875 


3.36124192571 


3.36124192571 


7 


-0.002317 


1.578558 


-0.00000693015 


3.36123499556 


15 


-0.000660 


1.577898 


-0.00000079070 


3.36123420487 


30 


-0.000183 


1.577715 


-0.00000009848 


3.36123410639 


60 


-0.000049 


1.577666 


-0.00000001246 


3.36123409393 


120 


-0.000010 


1.577656 


-0.00000000122 


3.36123409270 


extrap. 


-0.000015(15) 


1.577641(15) 


-0.00000000041(44) 


3.36123409229(44) 



in Sec. IHI CI it involves the derivative of the free elec- 
tron propagator with respect to the reference-state en- 
ergy. The reference-state contribution to the reducible 
part has been treated in Sec. IIIIB1 together with the 
reference-state contribution to the vertex term, thereby 
mutually cancelling the IR divergence inherent to both 
reference-state contributions. The total reference-state 
contribution is summarized in Eq. (|23|) . Left is the many- 
potential reducible part, 



A^ +) = (a\SV\a) 

x(a| 7 °^(E( £ )-E(°)(e)-S( a )( £ )) 



(54) 



Here, S^(e) and E(°)(e) are obtained from Eq. ([6]) by a 
replacement of the full Dirac-Coulomb Green function G 
by the free Green function and by the reference-state 
part of the propagator G^ a \ For the term with G^ a \ we 
have 

/oc 
duj (55) 
-OC 

xG (a »(e - w, x x ,x 2 ) a u D^(u, x 12 ) . 
In coordinate space, a representation of G^ a ' reads 

G y >(e - u;,x 1 ,x 2 ) = > — - , 56) 

z — ' £ — U — E a + lQ 



where we take into account all states with the same en- 
ergy as the reference state, i.e., also the state with op- 
posite parity but the same total angular momentum as 
compared to the reference state (pairs of states with the 
same |«| are energetically degenerate according to Dirac 
theory) . 

The evaluation of Eq. proceeds along the inte- 
gration contour Cm of P. J. Mohr [3(| for the (complex 



rather than real) photon energy. It is divided into a low- 
energy and a high-energy part. The low-energy contour 
C' L comprises the interval u> € (e a — i0, —iO) below the cut 
of the photon propagator and the interval (iO, e a + iO) on 
the upper bank of the cut, with e a being the reference- 
state energy. The high-energy contour C' H again consists 
of two intervals, [e a + iO, e a + ioo) and (e a — iO, e a — ioo). 
Because the low-energy part extends to comparatively 
high values of \oj\, the radial integrand for each single 
value of k become highly oscillatory. The behaviour of 
the integrand can only be improved if the full sum over in- 
termediate angular momenta is carried out before the ra- 
dial integrations. This is already evident from the model 
example given in Eq. (7.3) of Ref. (3l| . 



exp(-r[l - p\) 
r[l - p] 



- ^(2|«|+l)j w (ipr)/ig(ir), 
|k|=0 

(57) 

where j is a Bessel function and hS ' is a Hankel function 
of the first kind (0 < p < 1). The right-hand side of 
Eq. (|5"T|) involves functions that are highly oscillatory as 
a function of the radial variable r, but the left-hand side 
is a simple exponential. This "smoothing" phenomenon 
after the summation over the intermediate angular mo- 
menta is crucial for the evaluation as it enhances the 
rate of convergence of the multi-dimensional SE integrals 
dramatically. The convergence of the sum over |k| can 
be further accelerated by the so-called CNC transforma- 
tion [3l|. With maximum values of k in excess of 10 6 
being handled at ease using the CNC transformation, we 
are able to control the accuracy of the final evaluations. 
The derivative of the Green function is calculated directly 
using fourth-point and (alternatively, for verification) six- 
point difference schemes. We choose suitable values of 
the parameters so that the Green function derivative is 
calculated to a relative accuracy of 10~ 24 . Additional 
modifications are necessary in the extreme infrared re- 
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gion of photon energies; here the difference scheme is ad- 
justed so that the boundaries of the integration region are 
not crossed and sufficient accuracy is retained. A numer- 
ical subtraction of all singular terms due to lower-lying 
atomic states (e.g., the ground state) before doing any 
integrations over the photon energies and before evalu- 
ating the derivative of the Dirac propagator eliminates a 
potential further source of numerical loss of significance 
for the many-potential reducible part. 

G. Irreducible part 

With reference to Eq. ([5]) , we recall that the irreducible 
part is given as 

AE n = (H7°£(£a)|a) + (a| 7 °E(e a )|5a) , (58) 

with the renormalized SE operator E and the perturbed 
wave function [see Eq. Q] 



TABLE II: Individual contributions to the SE correction to 
the hfs of the IS state of hydrogen, in units AE/[a/ir Ep]. 
The specific contributions are discussed in Sec. IIII 51 (AE^ a '), 
Sec. ETC] (AE^J d and AE ( v °l), Sec. |mD] (AE ( V H), Sec. UJH] 
(A£<e r +) ), Sec.HmtAiw^), and Sec. |mG] (AE iT ). 



AE iT -0.01096549784 (5) 
AEl°l 8.28956864683 

AB re 1 d" > -3.83854412893 (5) 

AM2 -5.57958625925 

AEill -1.7835813412 (16) 
AEilt } 3.3612340923 (4) 

AE (a) -0.00002366906 
Total 0.4381018429 (16) 



\Sa)= Y, 

n 



i)(n\5V\c 



-a c n 



(59) 



We only need the diagonal-in-re component of the per- 
turbed wave function, because the SE operator is also 
diagonal in the total angular momentum. 

The evaluation of the irreducible part is carried out 
along the same contour Cm that is used for the many- 
potential reducible part. Within the high-energy part, 
the Green function is divided into two parts. The first 
is a subtraction term which involves a free propagator 
and an approximate one-potential term [3(| , which is ob- 
tained from the full one-potential term by commuting the 
Coulomb potential to the left of the electron propagators. 
The second is the remainder term which is the difference 
of the full and the approximate propagator. The subtrac- 
tion term contains all UV divergences of the irreducible 
part; these are cancelled against the mass counter term 
dm. The subtraction term is evaluated in momentum 
space, in a noncovariant integration scheme adjusted for 
bound-state calculations, where the spatial components 
of the photon momentum are integrated out before the 
photon energy integration. 



IV. NUMERICAL RESULTS 

Our calculation of the SE correction to the hfs and the 
g factor of hydrogen-like ions was performed in the Feyn- 
man gauge and for a point nucleus. The fine structure 
constant of a -1 = 137.036 was used in the calculation. 
The small deviation of this value from the currently ac- 
cepted one does not influence the numerical results 
for the higher-order remainder. An example set of indi- 
vidual contributions to the SE correction to the IS hfs 
of atomic hydrogen is presented in Table [Til 



o 
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FIG. 2: The higher-order remainder F„s{Za) for the SE cor- 



The SE correction to the hfs of an nS state can be 
represented as 



AE nS = -E F {nS) 



aoo + (Za) a w 



+ {ZaY{ \a*[(Za)-*\ a 22 + \n[{Za)- z ] a 2l + a 20 



(Za) 3 ln[(Za)- 2 } a 31 + (Zaf F nS (Za) 



(60) 



where Ep(nS) is the non-relativistic hfs value, and the 
Ojj are coefficients of the Za-expansion with the first 
index corresponding to the power of Za and the second 
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TABLE III: SE correction to the hfs of nS states of hydrogen- 
like ions. 5E n s = AE n s /[a/n Ep(nS)] and F n $ is the higher- 
order remainder defined by Eq. (|60[) . 



IS 


Z 


8E 1S 


F ls {Za) 




1 


0.4381018429 (16) 


-13.8308 (43) 




2 


0.373467600 (3) 


-14.1170 (9) 




3 


0.307583838 (4) 


-14.4121 (3) 




1 


0.241005731 (5) 


-14.6962 (2) 







U. 1 f ^ / J 


1 A Qfi7Q (0\ 

— 14. yO ( o \£) 






D 10681 5805 h 1 ^ 


— 15 9964 ('] ) 




7 

> 


03Q4764Q f 91 


— 1 5 4759 h 1 




g 


—0 097Q3933 (9) 


— 1 5 71 56 h 1 

J. O . I -L OU I -L 1 




g 


—0 09537946 (31 


— 1 5 9487 ("1 1 




i n 


— 1 698^*359 (4^ 


— 16 1 769 M ^ 

1U. J- 1 UZr ^ _1_ ^ 




11 


—0 2' : tfn57SQ 


—16 sooo (i ^ 

lu.OJi/W 1 11 




12 


— 9Q7Q0470 (A\ 


— 1 6 61 81 H S 


2S 


Z 


SE 2S 


F 2S {Za) 




1 


0.438692275 (3) 


-6.1205 (85) 




2 


0.375352042 (3) 


-6.9126(11) 




3 


0.311203194(5) 


-7.5833 (5) 




4 


0.246665425 (7) 


-8.1697(3) 




r. 



U.loiyooDoi [t)) 


— o. ( uoy i^z j 




u 


1 1 71 9^*309 (~\ 'W 


y.iy i o ^z. ^ 




7 


0^99646^ (9^ 


— 6^46 fl 




o 


— 01 9f>9 c iQ4 (91 


— in ns^Q h \ 




q 




— 10 4Q64 H ^ 




10 


—0 1 4955836 f3l 


— 10 8Q01 n 




11 


— 907661 5Q faS 


— 1 1 9701 n ^ 




12 


—0 979Q1 9S8 f4^ 


— 1 1 6SQ0 n 


3S 


Z 


SE 3S 


F 3S (^a) 




1 


0.43893143 (3) 


-1.727(69) 




2 


0.37613687(3) 


-2.6138 (90) 




3 


0.31274865 (3) 


-3.3530 (27) 






n 9ztQi 4i rc; (iS 


QQQzl (1 91 




5 


0.18548734(3) 


-4.5806 (7) 




6 


0.12186545 (3) 


-5.1133 (4) 




7 


0.05830586 (3) 


-5.6092 (2) 




8 


-0.00519144(3) 


-6.0757(1) 




9 


-0.06864589 (3) 


-6.5189(1) 




10 


-0.13209008 (3) 


-6.9430(1) 




11 


-0.19556594 (4) 


-7.3515(1) 




12 


-0.25912227(5) 


-7.7473(1) 



TABLE IV: SE correction to the hfs of 2Pj states of hydrogen- 
like ions. 5E n pj — AE n pj /[a/n Ep(nPj)] and G„p 7 is the 
higher-order remainder defined by Eq. (164[) . 



2Pl/2 


Z 


SE2P, 


G 2 p 1/2 (Za) 




1 


94939701 8 f5l 


— 1 1 393391 ^861 

_L_L. OZOOZ± IOUI 






9487 ("5l a 






2 


94801 6543 f5l 


—9 31 1 768 f95l 




3 


9460871 70 f7l 


—8 1 64978 ( 1 41 




4 


94371 9931 f7l 


— 7 370786 (Q) 




5 


0.240985405 (8) 


-6.771355 (6) 






0.2397 a 






6 


0.237932761 (8) 


-6.294696 (4) 




7 


0.234597810 (8) 


-5.902768 (3) 




8 


0.231007222 (8) 


-5.572857 (2) 




9 


0.227180996 (8) 


-5.290309 (2) 




10 


0.223134035 (8) 


-5.045123 (2) 






0.2202 a 






11 


0.218877214 (8) 


-4.830170 (1) 




12 


0.214418110 (9) 


—4.640191 (1) 


2P3/2 


Z 


5E 2 p 3/2 


G 2 p 3/2 (Za) 




1 


-0.12499329 (1) 


0.12609(18) 






-0.1254 6 






2 


-0.12498309 (1) 


0.079405 (55) 




3 


-0.12498458 (2) 


0.032176 (39) 




4 


-0.12501321 (2) 


-0.015499 (29) 




5 


-0.12508457(3) 


-0.063528 (21) 






-0.1255 6 






6 


-0.12521458 (3) 


-0.111933(16) 




7 


-0.12541922 (3) 


-0.160663(12) 




8 


-0.12571467 (3) 


-0.209698 (9) 




9 


-0.12611728 (3) 


-0.259028 (7) 




10 


-0.12664357(3) 


-0.308643 (5) 






-0.1271 6 






11 


-0.12731021 (3) 


-0.358538 (5) 




12 


-0.12813408 (3) 


-0.408711 (4) 



a Ref. hj 
6 Ref. [l|| 



020(15) = 17.12233875, 
020(25) = 11.90110542, 



a 20 (35) = 10.41704775, 



(63) 



corresponding to the power of the logarithm. We have 

o 00 (ri5) = 1/2 , oio(nS) = -8.03259003 , 

a 22 (nS) = - 2/3 , o 3 i(n5) = -13.30741592 , (61) 

o 2 i(15) = - 1.334503593, 

021(25) = 0.317103926, 

a 21 (35) = 0.921048823, (62) 



see the recent articles [32J [3J, |35j and references therein 
for earlier studies. F n s is the higher-order remainder, 
which we address in our numerical all-order approach. 
Our numerical results for the SE correction to the hfs of 
the 15, 25, and 35 states are listed in Table ITTTT 

The SE correction to the hfs of an nPj states is much 
less studied. Only the leading term of its Za expansion 
is known today. The correction, therefore, is written as 



AE nPj =E F (nPj) 



a a + (Za) 2 G nPj {Za) , (64) 



with G n p, being the higher-order remainder. The coeffi- 
cient Qqo is given by aoo^-Pi/^) = 1/4 and 000(^-^3/2) = 
— 1/8 [36]. Our numerical results for the SE correction 



12 




-12- 



— i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i— 
01 23456789 10 11 12 

z 

FIG. 3: The higher-order remainder G n pj(Za) for the SE 
correction to the hfs of the 2Pi/ 2 and 2P 3 / 2 states. 
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FIG. 4: The higher-order remainder for the SE correction to 
the g factor of the 2P 1/ / 2 and 2P 3 / 2 states. 



TABLE V: SE 


correction to the g-factor of nS states of 


hydrogen-like ions, in ppm. H n s is the hij 


'her-order remain- 


der. 






lb Z 




His(Za) 


1 


2322.840230 (2) 


19.(3.) 


2 


2322.904037 (4) 


23.3 (2.6) 


3 


2323.014295 (8) 


22.88 (70) 


4 


2323.175525 (14) 


22.58 (29) 


r 




OQOO QflOflQ t o\ 

zozo.o9/yo \Z J 


22.36 (15) 


6 


2323.67243 (3) 


OO 1 Q [Q\ 


7 


2324.02001 (5) 


22.02 (6) 


8 


2324.44213 (7) 


21.87(4) 


9 


2324.94538 (8) 


21.71 (3) 


10 


2325.53651 (10) 


21.57(2) 


11 


2326.22235 (13) 


21.42 (2) 


12 


2327.00983 (12) 


21.28(1) 


O C V 


&92S 




1 


2322.824624 (2) 




2 


2322.840323 (8) 




3 


2322.86696 (2) 


19. (11.) 


4 


2322.90509 (3) 


22.1 (4.7) 


r 

O 


OQOO CiKKQQ f K\ 

/ozz.yoooo (o ) 


22.5 (2.4) 


6 


2323.01836 (6) 


OO z. (~\ 1\ 
ZZ.O [L.o ) 


7 


2323.09490 (8) 


22.43 (74) 


8 


2323.18571 (8) 


22.31 (42) 


9 


2323.29154 (9) 


22.18 (24) 


10 


2323.41319 (9) 


22.05 (15) 


11 


2323.55142 (11) 


21.92(11) 


12 


2323.70704 (12) 


21.79 (8) 


3S Z 




H 3S (Za) 


1 


2322.821746 (7) 




2 


2322.828684 (10) 




3 


2322.84038 (2) 




4 


2322.85698 (3) 


19. (18.) 


5 


2322.87866 (5) 


20.7(8.4) 


6 


2322.90560 (6) 


21.5(4.6) 


7 


2322.93798 (9) 


21.8(2.9) 


8 


2322.97600 (12) 


22.0(2.0) 


9 


2323.0198 (2) 


22.0(1.4) 


10 


2323.0697 (2) 


22.0(1.0) 


11 


2323.1258 (2) 


21.89(77) 


12 


2323.1882 (3) 


21.81(59) 



to the hfs of the 2Pi/ 2 and 2P 3 / 2 states are listed in Ta- 
blelVl 

The results for the higher-order remainders F n s and 
G n p r inferred from our numerical data are plotted in 
Figs. [2] and [31 respectively. For the 2P 1 / 2 state, a fit of 
our results is consistent with the Za expansion of the 
form 

G n p 1/2 (Za) = a 2 i ln[(Zay 2 ] + a 2 o + ■ ■ ■ , (65) 

where the value of the logarithmic coefficient is very close 
to a2i(2Pif 2 ) = —3/2 and the constant term is about 
a 20 (2P 1 / 2 ) = 3.5. For the 2P 3 / 2 state, the numerical 
data are consistent with a 2 i(2P 3/ / 2 ) = 0. Our results in 
Table ITVl are in moderate agreement with those obtained 



previously in Refs. (l3l.[l3| but significantly improve upon 
them in numerical accuracy. Nevertheless, we disagree 
with the suggestion [l3| about the possible presence of 
the squared logarithm in the Za expansion (|65[) for the 
2P x / 2 state. A more careful investigation of the analytic 
structure of the hig her-order terms is performed in the 
follow-up paper [37] . 



The SE correction to the bound-electron g factor of an 
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TABLE VI: SE correction to the g factor of 2Pj states of 


hydrogen-like ions, 


in ppm. I n pj is the higher 


-order remain- 


der. 






op , 7 


A 52P 1/2 




1 


-774.258151 (3) 


0.121258 (21) 


2 


-774.212929 (11) 


0.121715 (22) 


3 


-774.13687(2) 


0.122414(19) 


4 


-774.02917(3) 


0.123280 (14) 


5 


-773.88876 (3) 


0.124305 (10) 





77Q 71 /I/19 /''^ 


U. 1Z04 i ) 


7 


-773.50460 (3) 


0.126803 (5) 


8 


-773.25838 (3) 


0.128186 (4) 


9 


-772.97356 (3) 


0.129711 (2) 


10 


-772.64862 (3) 


0.131336 (2) 


11 


-772.28175 (4) 


0.133054(3) 


12 


-771.87108 (8) 


0.134858 (5) 


op„,„ 7 


A 52P 3/2 




1 


774.291470 (3) 


0.148104 (21) 


2 


774.346522 (12) 


0.148294(24) 


3 


774.43854 (3) 


0.148567(24) 


4 


774.56774 (4) 


0.148851 (22) 


5 


774.73495 (6) 


0.149338 (18) 


6 


774.94028 (6) 


0.149816(14) 


7 


775.18442 (6) 


0.150350(11) 


8 


775.46799 (6) 


0.150933(7) 


9 


775.79167(5) 


0.151561(4) 


10 


776.15615 (5) 


0.152231 (4) 


11 


776.56218 (6) 


0.152940 (4) 


12 


777.01055 (6) 


0.153685 (4) 



nS state can be represented as 
(Za) 2 



A a 

Ag n s = - 

7T 



1 



-I^(ln[(Za)- 2 ]6 41 +6 40 
^— U nS {Za) 



(66) 



where the bij are known coefficients of the Za expansion: 



b 2 o(nS) = - , &41 (nS) = 
6 

640(15) = - 10.23652432, 
b i0 (2S) = - 10.70771560, 
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640(35) 



11.52963397, 



(67) 



see Ref. [38| and references therein. H n s is the remain- 
der incorporating all higher-order contributions. It is re- 
markable that the higher-order remainder H n s enters in 
the relative order (Za) 5 rather than in the relative or- 
der (Za) 3 , as in the case of the hfs. This means that 
cancellations in extracting the remainder from numerical 
results for Z = 1 are larger for the g factor than for the 
hfs by four orders of magnitude. 



Our numerical results for the SE correction to the g 
factor of the electron 15, 25, and 35 states of light 
hydrogen-like ions are presented in Table |V] We observe 
that the higher-order remainder behaves very similarly 
for all nS states studied, the 25 remainder being just 
about 2% larger than that for the 15 states and the 35 
and 25 remainders being equal within the numerical un- 
certainty. The accuracy of the direct numerical determi- 
nation of the 15 remainder for Z = 1 and Z = 2 can 
easily be increased by extrapolating values obtained for 
higher values of Z . An extrapolation yields the improved 
results His (la) = 23.39(80) and H is (2a) = 23.03(44). 
Improved values of the higher-order remainder for the 
25 and 35 states are most easily obtained by scaling the 
15 remainder. The trend of the higher-order remain- 
der for low Z is consistent with a numerically large, n- 
independent coefficient b 50 in Eq. (|66|) . 

For the Pj states, the bound-electron g factor is stud- 
ied less thoroughly than for the 5 states. The leading 
term of its Za expansion is due to the electron anomalous 
magnetic moment (amm) and is immediately obtained for 
a general state as [39| 

1 - 2k 

^WTT)' (68) 

where k is the Dirac quantum number and j is the to- 
tal angular momentum of the electron state. For the P 
states, the explicit results are 600(^^1/2) = ~ 1/3 and 
boo(nP 3/2 ) = 1/3. 

The next-order term, 620 , consists of two parts, one in- 
duced by the electron amm and the other, by the emission 
and the absorption of virtual photons of low energy (com- 
mensurate with the electron bin ding energy). A simple 
calculation of the first part gives [3!| 620(^^1/2, amm) = 
— 1/2 and 620(^^3/2; amm) = 1/10. The second part is 
nonvanishing for states with I 7^ only and is more com- 
plicated. Its general expression is known Hl| but 
the only numerical result available for hydrogenic atoms 
is the estimate made in Ref. [42j], which disagrees with 
our numerical values both in the sign and the magni- 
tude. Commenting on this fact, we note that the esti- 
mate is based on a rather crude approximation. Namely, 
the sum over the entire discrete and continuous spectrum 
of virtual states was replaced by the contribution of the 
lowest 12 discrete bound states only. We argue that such 
approximation might be inapplicable for the problem in 
hand. The reason is that, e.g., for the Bethe logarithm 
(which is also a contribution induced by the low-energy 
photons) the dominant contribution originates from the 
continuum spectrum [43[ so that such approximation is 
clearly inadequate. 

Since the (Za) 2 term is not presently known anyalyt- 
ically, we define the higher-order remainder for the Pj 
states as 



(Za) 2 I nPj (Za) 



(69) 



Our numerical results for the SE correction to the g 
factor of the electron 2^/2 and 2P 3 / 2 states of light 
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hydrogen-like ions are listed in Table IVII The corre- 
sponding higher-order remainder function is plotted in 
Fig. m 



V. CONCLUSIONS 

We have discussed, in detail, a numerical evaluation, 
nonperturbative in the binding Coulomb field, of the self- 
energy correction to the hyperfine splitting and of the 
self-energy correction to the g factor in hydrogen-like 
ions with low nuclear charge number Z = X, ... ,12. We 
consider the ground state, as well as the 2S and 35 ex- 
cited states, and the 2P 1 / 2 as well as 2P 3 / 2 states. The 
value of a -1 = 137.036 is employed in all calculations. 
At the level of precision we are operating at, the final 
results for the self-energy corrections depend on the pre- 
cise value employed very sensitively. However, the main 
dependence on the value of a is accounted for by the an- 
alytically known lower-order terms. Thus, the results for 
the remainder functions F n s(Za), G n Pj(Za), H„s(Za), 
and I n pj (Za) as given in Tables InTlUVIIVl and lVII are not 
influenced by the value of a employed. Even if a value 
of a which differs from aT 1 — 137.036 on the level 10 -7 
were employed, then the values of the remainder func- 
tions would not change: their main uncertainty is due to 
limits of convergence of the integrals that constitute the 
nonperturbative self-energy corrections, as described in 
the preceding sections of this article. 

The organization of our calculation is described in 
Sec. IIIIl We consider separately the reference-state con- 
tribution to the reducible and the vertex part (Sec. lIIIB]) 
and the zero-potential contribution to the vertex and 
to the reducible part (Sec. IIIIC|) . The one-potential 
and many-potential vertex parts, which represent the 
most challenging part of the calculation, are discussed in 
Sees. MI Dl and [HI El The many-potential reducible part 
and the irreducible part conclude the discussion of our 
computational method fSecs. lHTFl and lHI G[) . Numerical 
calculations were carried out on the parallel computing 
environments of MPI Heidelberg and MST Rolla. 

It is instructive to compare the numerical results ob- 
tained (Tables IIIIl — IVI|) to analytic results from the Za 
expansion. The analytic parameterization of the self- 
energy correction to the hyperfine splitting according to 
Eq. (foT))) entails both logarithmic as well as nonlogarith- 
mic corrections. Our numerical results for the scaled 
self-energy correction SE n to the hyperfine splitting and 
for the nonperturbative remainder function F n s(Za) are 
given in Table IIIIl for S states. A plot of the data 
(see Fig. ^ indicates that the higher-order remainders 
F n s(Za), for Z —* 0, may converge toward an 030 coef- 
ficient which is significantly dependent on the principal 
quantum number. 

The scaled self-energy correction 5E n pj for 2Pi/ 2 and 
2P 3 / 2 states is analyzed in Table IIVI A plot of the 
data (see Fig. aids in a comparison to an analytic 
model for the correction, given in Eq. (|65p . The non- 



perturbative remainder G n pj {Za) for the hfs of P states 
is seen to be well represented by an analytic model of 
the form {021 ln[(Za)~ 2 ] + a 2 Q + ...}, where a fit to the 
numerical data indicates that a2i(2-Pi/2) = —3/2 and 
that a2i(2P 3 / 2 ) = 0. It has been suggested in Ref. [l3| 
that a "double" (squared) logarithm in the Za expan- 
sion could be present for low nuclear charge; the latter 
would correspond to a nonvanishing 022 coefficient for P 
states. Our numerical results in Table HVl do not contra- 
dict those of Ref. [l3[ on the level of numerical accuracy 
obtained in the cited reference. However, we cannot con- 
firm the presence of such a double-logarithmic correction 
(see also [37| )• 

A large number of analytic terms are known for the 
self-energy correction to the g factor for S states [see 
Eq. (|66p ]. The higher-order remainder H„s(Za) for the 
g factor of S states is thus "separated" from the leading- 
order effect by about ten orders of magnitude for Z = 1. 
Thus, although our direct numerical evaluation of the 
self-energy correction for the ground state at Z = 1 is 
precise [Ag is = 2322.840230(2) x 10~ 6 at Z = 1], we 
can only infer the higher-order remainder H\s(Za) at 
Z = 1 to about ±10%: the result after the subtraction of 
lower-order terms is His (la) = 19 (3). By extrapolation 
of more accurate data for the remainder obtained from 
higher values of the nuclear charge, we can obtain the 
improved results Hi S (la) = 23.39(80) and H is (2a) = 
23.03 (44). The remainder functions at very low Z appear 
to depend only very slightly on the principal quantum 
number; they are consistent with H n s(Za) approaching 
an n-independent coefficient &50 as Z — > 0. 

For the g factor of P states, only the leading coefficient 
&00 is known from the Za expansion [see Eq. (|69p ]. The 
self-energy remainder function I n Pj (Za) for the g factor 
of P states can be inferred from our numerical data in 
Table IVII after subtraction of the leading analytic term 
as given in Eq. (|69p. The numerical data for P states are 
consistent with the functions I n p r (Za) tending toward a 
constant for Z — > 0. A plot of the data in Fig. 2] confirms 
this trend. 

To conclude, we have performed an all-order (in Za) 
calculation of the self-energy correction to hyperfine split- 
ting and g factor in hydrogen-like ions with low nuclear 
charge numbers. The calculation is accurate enough 
to infer higher-order remainder terms without any ad- 
ditional extrapolation, by a simple subtraction of the 
known terms in the .Za-expansion. We improve the nu- 
merical accuracy by several orders of magnitude as com- 
pared to previous evaluations; this leads to improved the- 
oretical predictions for all QED effects considered in this 
article. 
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that X is a combination of the Legendre polynomials 
with some coefficients. It is straightforward to find the 
coefficients by performing integrations in Eq. (|A4|) ana- 
lytically for each particular case (where advantage may 
be taken of computer algebra). 



APPENDIX A: ANGULAR INTEGRATIONS IN 
MOMENTUM SPACE 

In this section we demonstrate how to perform the in- 
tegration over the angular variables in momentum space. 
The problem in hand can be formulated as follows: the 
general expression of the form 



J dpidp 2 F(p lr ,p 2r ,^)A(pi)B(p 2 ) 



(Al) 



where F, A, and B are some arbitrary functions and £ = 
pi ■ pi , needs to be integrated over all angular variables 
exept for £, i.e., to be reduced to the form 

J d£F{p lr ,p 2r ,£)X{A,B-t). (A2) 

In order to write a general expression for the function 
X(£) in terms of A and B, we use the standard decom- 
position of the function F in terms of the spherical har- 
monics Yi m , 

F(pir,P2r,0 = 2irJ2YL(Pl)Ylm(P2) 
I III 

x J d^F{ Plr ,p 2r ,0Pi(0> (A3) 

where Pi are the Legendre polynomials. From this, we 
immedately have 



OO I r „ 

X(A,B;0 =2ttX>(0 E / d P* A (Pi) Y L(Pi) 



1=0 m=-l 
dp 2 B(p 2 )Yl m (p2 



(A4) 



This general formula greatly simplifies when one of the 
functions (say, B) is unity or the identity function, i.e., 
B{p 2 ) — p 2 . When B is unity, only the term with I = 
contributes, and we have: 

X(A,l;£) = 2n[dpA(p). (A5) 



When B = id, 



X(A,id;0 = 2< J dppA(p). 



(A6) 



In more complicated cases with B — pipk ■ ■ ., formulas 
for X can be in priciple obtained by using the Racah 
algebra. Alternatively, one can observe [from Eq. (|A4[) ] 



In the present work, we need angular integrals of three 
types, Ki, K 2l and K 3 , defined as (^ = 1/2): 



3i f 

— / dp x dp 2 F(p lr ,p 2r ,(,)xl^(Pi) [Pi x <r] z X- K AP2) 

= J F(pi r , par, O-^i («). ( A7a ) 

— / dpxdp 2 F(p lr ,p 2r ,(,)xi, li (pi) [P2 x (t] z X- k AP2) 

= J ^F{p lriP2ri 0K[(K), (A7b) 

3i f 

— / dp x dp 2 F(p lr ,p 2r ,(,)xl }l ,(pi) [pi XftUvfe) 

= J dfiF(p lr ,p 2r) 0K 2 (K), (A7c) 

— / dpidp 2 F( Plr ,p 2r ,£,)xl^(Pi)i(TzX K AP2) 

= J ^F(p lr ,p 2r ,Z)K 3 (K). (A7d) 



Using the technique described above, we obtain the fol- 
lowing results for the basic angular integrals: 



Ki(k) 



K' 1 {k) = 



K = -l, 
1 , K = 1 , 

|(l-3£ 2 ), « = -2, 

U> K = 2 > 

-1, K=-l, 
L K = 1 , 

-§£, k = -2, 
[ -I(1-3C 2 ), « = 2, 



(A8a) 



(A8b) 
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as well as 



K 2 (k) 



-i(i-a, «=-2, 



•&£(i-£ 2 )> « = 2, 



(A8c) 



if 3 («) 



_ 3 
2 ' 



K = — 1 . 
K = 1 , 
K= -2. 
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(A8d) 
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